## calculate Social Cost of Carbon
## - damages from DICE models
## - term structures of SDRs from different time series models

rm(list=ls())
graphics.off()
library(dplyr)

##################################################
## CHOOSE HERE

## Table 2, top panel
## damages from Haensel et al DICE model
cat("# Haensel-et-al DICE model\n")
d <- read.csv("data/dice_haensel_damages.csv") %>%
    filter(year <= 2410)
mats <- d$year - d$year[1]
loss <- d$damages
cat("SCC model:", crossprod(d$damages, d$DF), "\n")

## ## Table 2, top panel
## ## damages from DICE-FAIR-Geoffroy model in Dietz et al
## cat("# DICE-FAIR-Geoffroy model from Dietz et al\n")
## d <- read.csv("data/dice_dfg_damages.csv") %>%
##     filter(year <= 2410) # max 400 years
## mats <- d$year - d$year[1]
## loss <- d$damages
## cat("DFG model (published 17.45):", crossprod(d$damages, d$DF), "\n")

## APPENDIX RESULTS

## ## damages from DICE-96 model -- via Newell-Pizer
## cat("# DICE-1996 model\n")
## d <- read.csv("data/dice_newell_pizer_damages.csv")
## loss <- d$loss # marginal damages
## mats <- 0:400 # maturities/horizon of future damages

## ## damages from DICE-2016 model
## cat("# DICE-2016 model\n")
## ## - 5 years per period, but damages account for that (no need to scale)
## d <- read.csv("data/dice_2016_damages.csv") %>%
##     filter(year <= 2410)  # max 400 years
## mats <- d$year - d$year[1]
## loss <- d$damages
## cat("Nordhaus (PNAS published value 31.2):", crossprod(d$damages, d$DF), "\n")

##################################################

stopifnot(length(loss) == length(mats))
cat("SCC with 4-percent discount rate:", round(crossprod(1.04^(-mats), loss),1), "\n")
cat("SCC with 2-percent discount rate:", round(crossprod(1.02^(-mats), loss),1), "\n")

## six different time series models
files <- c(
           "results/sdr_uc_y1.RData",
           "results/sdr_uc_y10.RData",
           "results/sdr_meanshift_y1.RData",
           "results/sdr_meanshift_y10.RData",
           "results/sdr_ar3_expon_y1.RData",
           "results/sdr_ar3_expon_y10.RData")
J <- length(files)
nms <- c("rstar change", "SCC-1990", "SCC-2020", "rel. change")
tbl <- matrix(NA, J, length(nms))
## names of time series models
rownames(tbl) <- c(
    "UC model, 1y rate",
    "UC model, 10y rate",
    "AR model, break, 1y rate",
    "AR model, break, 10y rate",
    "AR model, learning, 1y rate",
    "AR model, learning, 10y rate")
colnames(tbl) <- nms
for (j in 1:J) {
    load(files[j])
    N <- length(P1)
    stopifnot(N == 400)
    P1 <- c(1, P1)
    P2 <- c(1, P2)
    SCC1 <- crossprod(P1[mats + 1], loss)
    SCC2 <- crossprod(P2[mats + 1], loss)
    tbl[j,] <- c(rstar_new-rstar_old,
                 SCC1,
                 SCC2,
                 100*(SCC2/SCC1-1))
}

print(round(tbl,1))
